% 9-March-2007 Note: Chnages have been made to the orginal file
% eej_induction_paper.m
% The numerics in the paper EEJinduction.doc belong to the earlier file

% Process the eej strength derived from satellite data
% and compare it with ground deribed eej indices
%start date 30 OCT 2005
%latest date 23 SEP 2006
% The idea is just to see the slop of the regression at 0-degree bin
%clear;
more off;
load c:\manoj\projects\eej\OBS_array.mat OBS
load c:\manoj\projects\eej\mjdkpf107 kp;

%OBSdH1 description 
%The file OBS is Structured array with stations as listed below
%[1-ABG,2-ETT,3-HYB,4-TIR,5-HUA,6-FUQ,7-MBO,8-GUI,9-AAE,10-QSB,11-GUA,
%12-CBI,13-PND,14-ELT,15-AAE(new),16-NGP,17-UJJ,18-KOD,19_YAP,20-BIK]

% Note please update c:\manoj\projects\eej\used_obs.xls if more data are
% added to OBS_array.mat

%-------- correlate the CHAMP derived EEJ strength and ground delta H
%--------------------------------------------------------------------
fidd = fopen('c:\manoj\temp\eej_slope_new.txt','at');

%obs_pairs = [2,3;4,1;5,6;19,20;15,14];
obs_pairs = [2,3;4,1;5,6;19,20;9,14]; %Olde AAE data - ELT
%Note NGP !
tic;
%Control Center
%[1-ABG,2-ETT,3-HYB,4-TIR,5-HUA,6-FUQ,7-MBO,8-GUI,9-AAE,10-QSB,11-GUA,
%12-CBI,13-PND,14-ELT,15-AAE(new),16-NGP,17-UJJ,18-KOD,19-YAP, 20-BIK ]

eej_data_set    = 3;%1 -CHAMP Scalar Version 1, 
                    %2 -CHAMP Scalar Version 2
                    %3 -CHAMP SCalar Version 3, 
                    %4 -CHAMP vector        
                    %5  -Oersted scalar
                    %6  -Oersted vector
                    %7 - SAC-C scalar
                    %8 - All scalar
satdata  = 2; % %Pear Current Density = 1; av cur distan 2; %TWO PLEASE !!
data_set        = 3;% 1-common; 2-Old-independent; 3 old - all (The last two options should make similar results)

obslt    = [9.9 13.1]; %NOTE the LT are represented as for eg 8.99999 or 9.001 the filter
incrlat  = 10;

mindate  = tmjd(2000,01,01,0,0,0,1); % starting date
maxdate  = tmjd(2002,12,31,0,0,0,1)+1; % ending date

robft    =      1; %Robust outlier elimination and slope determination
Rejection_factor = 0.1;% 0.1 to 0.9 -> amout of data to be rejected (do not put this to 1 !)
plotcol  =      ['k','k','k','k','k','k'];
kp_thrashold =  2;
dummy_p =       1;
fprintf(fidd,'eej_data_set %d, data_set %d, satdata %d, robft = %d\n',eej_data_set, data_set, satdata, robft);
fprintf('eej_data_set %d, data_set %d, satdata %d, robft = %d\n',eej_data_set, data_set, satdata, robft);

switch eej_data_set,
    case 1,
        switch data_set,
            case 1,
                load c:\manoj\projects\eej-induction\common_passes common_old_data;
                data = common_old_data;
            case 2,
                load c:\manoj\projects\eej-induction\common_passes used_old_data;
                data = used_old_data;
            case 3,
                data = load('c:\manoj\projects\eej\eej_par_all.txt');
        end;
            switch satdata,
                case 1,
                eej_index = 9;
                case 2,
                eej_index = 10;
            end;
    case 2,
        switch data_set,
            case 1,
                load c:\manoj\projects\eej-induction\common_passes common_new_data;
%                 load c:\manoj\projects\eej-induction\bad_passes L;
%                 data = common_new_data(~L,:);
                  data = common_new_data;
            case 2,
                load c:\manoj\projects\eej-induction\common_passes used_new_data;
                data = used_new_data;
            case 3,
                load c:\manoj\projects\eej\EEJ_Stefan data;
            end;
        switch satdata,
            case 1,
            eej_index = 5;
            case 2,
            eej_index = 6;
       end;
    case 3,
        data = load('c:\manoj\projects\eej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\CHAMP_scalar.txt');
        switch satdata,
            case 1,
            eej_index = 5;
            case 2,
            eej_index = 6;
        end;
    case 4,
        data = load('c:\manoj\projects\eej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\CHAMP_vector.txt');
        switch satdata,
            case 1,
            eej_index = 5;
            case 2,
            eej_index = 6;
        end;

    case 5,
        data = load('c:\manoj\projects\eej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\Oersted_scalar.txt');
        switch satdata,
            case 1,
            eej_index = 5;
            case 2,
            eej_index = 6;
        end;
        
    case 6,
        data = load('c:\manoj\projects\eej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\Oersted_vector.txt');
        switch satdata,
            case 1,
            eej_index = 5;
            case 2,
            eej_index = 6;
        end;
        
    case 7,
        data = load('c:\manoj\projects\eej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\SAC-C_scalar.txt');
        switch satdata,
            case 1,
            eej_index = 5;
            case 2,
            eej_index = 6;
        end;
    case 8,
  data = load('c:\manoj\projects\eej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\SACOERSTEDCHAMP_Scalar.txt');
% [y,i] = sort(data(:,1));%This is necessary since the data file is just stitch together
% data = data(i,:);%of CHAMP, Oersted and SAC-C files - > not
% necessary ! 
% data(1:8718,:) = []; %Removing sac-c scalar data(SAC-C upto number 1:8718)
           
        switch satdata,
            case 1,
            eej_index = 5;
            case 2,
            eej_index = 6;
        end;
    
end;
%INFORMATION
%The EEJ parameters in the file
%1-mjd,2-longitude of eq-crossinfg, 3-zero of current position-south
%4 -north,5-half-width of current - south,6-north
%7minimum current south,8-north,9-pear current density 10-average
%difference 11 - diff to south minumum 12 - diff to north minum, 13-total
%east current ,14 total return current., 15 -f,16 - eq angle, 17-
%eq-direction, 18 - mgnteic equator latitude

for ikk = 1:length(obs_pairs),

    %     if ikk == 5,
%         robft = 0; %This is becasue for AAE-ELT the best slope obtaine without robust fitting !
%     end;
eq_obs = obs_pairs(ikk,1);
sq_obs = obs_pairs(ikk,2);
data1  = [];    
    

        obs_index = 2;
        data1(:,1) = OBS(eq_obs).mjd;
        data1(:,2) = OBS(eq_obs).data - OBS(sq_obs).data;
        obslon = OBS(eq_obs).long;
        if obslon > 180,
            obslon = -(360-obslon);
        end;
        st = obslon - incrlat/2;
        en = obslon + incrlat/2;


L = data(:,2) >= st & data(:,2) <= en & data(:,1) < maxdate  &  data(:,1) > mindate; %The last filter is to select only data before Dec 31, 2002
IndiaE = data(L,:);
ndata = 1;
ind1 = [];
ind2 = [];
nnodata=1;
for i = 1:length(IndiaE(:,1)),
   [trash ind] = min(abs(data1(:,1) - IndiaE(i,1))); %This is 2 times faster than findnearest 27sep06
  % ind = findnearest(IndiaE(i,1),data1(:,1));
    [lt1,dummy1] = champ_lt(IndiaE(i,1),IndiaE(i,2),obslon);
    [dummy2,lt2] = champ_lt(data1(ind,1),IndiaE(i,2),obslon); % This is because observatory FDAY ~= CHAMP FDAY
    
    if eej_data_set == 1,
            dummy_mjd(dummy_p) = IndiaE(i,1);
            dummy_p = dummy_p+1;
         if lt2 >= obslt(1) & lt2 <= obslt(2) & abs(data1(ind,obs_index)) < 300 & data1(ind,1) < max(data(:,1));% & kp(ind) <= kp_thrashold, %The last fliter take care the nodata (9999) points
             ind1(ndata) = ind;

             ndata = ndata+1;
        else,
            ind2(nnodata) = i;
            nnodata = nnodata + 1;
        end;
    elseif eej_data_set >= 2, %Format change
            dummy_mjd(dummy_p) = IndiaE(i,1);
            dummy_p = dummy_p+1;

            kp1 = IndiaE(i,10);
            cmax = IndiaE(i,5);
            if lt2 >= obslt(1) & lt2 <= obslt(2) & lt1 >= obslt(1) & lt1 <= obslt(2) ...
            & abs(data1(ind,obs_index)) < 300 & data1(ind,1) < max(data(:,1)) ...
            & kp1 <= kp_thrashold & cmax >= 0.03, %The second last fliter take care the nodata (9999) points

            ind1(ndata) = ind;
            ndata = ndata+1;
            else,
            ind2(nnodata) = i;
            nnodata = nnodata + 1;
            end;
      end;
end;

if length(ind2) > 0,
    IndiaE(ind2,:) = [];
end;

if length(ind1) > 1,
    obsE = data1(ind1,:);
    X = IndiaE(:,eej_index);
    Y =   obsE(:,obs_index);
    sel_mjd = IndiaE(:,1);
   if robft == 1&length(X)>3,    
%    Robust outlier rejection  -- START
     [rs,rstat] = robustfit(X,Y);
     residuals = abs(rstat.resid);
%      LLL = residuals == max(residuals); %Just remove the most noisy point
     [y,i] = sort(residuals);
     N_factor = round(length(residuals)*Rejection_factor);
     LLL = i(end-N_factor:end);
     
     X(LLL) = [];
     Y(LLL) = [];
     sel_mjd(LLL) = [];
     [rs,rstat] = robustfit(X,Y);
     c(1) = rs(2);
     c(2) = rs(1);
     % Robust outlier rejection  --- END
   else,
     [c,sfit] = polyfit(X,Y,1);
   end;
    [cof,err] = corrcoef(X,Y);  % correlation of data 

    %----------- PLOT SCATTER
    subplot(3,2,ikk);
    m = polyval(c,[min(X),max(X)]);
    plot(X,Y,[plotcol(ikk) 's'],'MarkerFaceColor',plotcol(ikk));
    hold on;
    set(gca,'FontSize',16);
    grid on
    h=xlabel('A/m');
    set(h,'FontSize',16);
    ylabel('\DeltaH');
    plot([min(X),max(X)],m,[plotcol(ikk) '-'],'LineWidth', 2);
    text(0,25,sprintf('%s-%s dH = %5.1f * I + %5.2f',OBS(obs_pairs(ikk,1)).code,OBS(obs_pairs(ikk,2)).code,...
        c(1),c(2)),'FontSize',16,'Color',plotcol(ikk));
    axis([0 0.3 0 200]);
    title(sprintf('%s-%s', OBS(obs_pairs(ikk,1)).code,OBS(obs_pairs(ikk,2)).code));
    ndd = length(X);

clear ind f_ch f_ob ind1 ind2 ndata nnodata IndiaE obsE ltt1 ltt2 ltt3 ltt4 obst chpt X Y;

end;

fprintf('\n');
fprintf('%s-%s %6.2f %6.2f %6.2f ', OBS(obs_pairs(ikk,1)).code,OBS(obs_pairs(ikk,2)).code,...
    c(1),c(2),OBS(obs_pairs(ikk,1)).eej_factor3*c(1));
fprintf(fidd,'%s-%s %6.2f %6.2f %6.2f ', OBS(obs_pairs(ikk,1)).code,OBS(obs_pairs(ikk,2)).code,...
    c(1),c(2),OBS(obs_pairs(ikk,1)).eej_factor3*c(1));

fprintf('%3d %4.2f %e ',ndd,cof(2,1),err(2,1));
fprintf(fidd,'%3d %4.2f %e ',ndd,cof(2,1),err(2,1));

fprintf(fidd,'\n');

end;
more on;
fclose all;
fprintf('\n');
toc;

%THE FOLLOWING SCRIPT IS USED TO PLOT THE CC AND SLOPE parameters
% Used in the figures c:\manoj\projects\eej-induction\CC*.fig and
% SLOPE*.fig

% legend1 = legend(gca,{'CH\_SC.1','CH\_SC.2','CH\_SC.3','CH\_VEC.3','OR\_SC.3','OR\_VEC.3','SAC-C\_SC.3','ALL\_SC.3'},'FontSize',12,'Position',[0.1513 0.1268 0.275 0.4048]);
% data_version = {'CH\_SC.1','CH\_SC.2','CH\_SC.3','CH\_VEC.3','OR\_SC.3','OR\_VEC.3','SAC-C\_SC.3','ALL\_SC.3'};
% 
% plot(a(:,1:2:end)','LineWidth',2)
% data_version = {'CH\_SC.1','CH\_SC.2','CH\_SC.3','CH\_VEC.3','OR\_SC.3','OR\_VEC.3','SAC-C\_SC.3','ALL\_SC.3'};
% legend([b(2:end,1)])
% axis([0 9 0 1000])
% set(gca,'XTickLabel',[])
% text([1:8],ones([1,8])*700,data_version,'rotation',90,'FontSize',8, 'FontWeight','bold')
% title('Slope of regression : \DeltaH Vs. ACD for KP <= 2, LT 10-13, LONG \pm 5^{o}')
% ylabel('Slope of regression nT.A^{-1}.m')
% % 
% plot(a(:,2:2:end)','LineWidth',2)
% axis([0 9 0 1.3])
% set(gca,'XTickLabel',[])
% legend([b(2:end,1)])
% data_version = {'CH\_SC.1','CH\_SC.2','CH\_SC.3','CH\_VEC.3','OR\_SC.3','OR\_VEC.3','SAC-C\_SC.3','ALL\_SC.3'};
% text([1:8],ones([1,8]),data_version,'rotation',90,'FontSize',8, 'FontWeight','bold')
% title('CC from \DeltaH Vs. PCD for KP <= 2, LT 10-13, LONG \pm 5^{o}')
% ylabel('CC')